Label transfer EDA

library(Seurat)
library(SnapATAC)
library(tidyverse)
library(DescTools)  # 4 AUC function

Registered S3 method overwritten by 'DescTools':
  method         from 
  reorder.factor gdata

library(glue)

Attaching package: ‘glue’

The following object is masked from ‘package:SummarizedExperiment’:

    trim

The following object is masked from ‘package:GenomicRanges’:

    trim

The following objects are masked from ‘package:IRanges’:

    collapse, trim

The following object is masked from ‘package:dplyr’:

    collapse

library(ggalluvial)  # 4 river plot
library(ggpubr)

model.cca <- readRDS("~/models/modelCCA_union_hvg_F74_SCElist_20191113.RDS")
model.liger <- readRDS("~/models/modelLiger_union_hvg_F74_SCElist_20191113.RDS")
model.conos <- readRDS("~/models/modelConos_union_hvg_F74_SCElist_20191113.RDS")
seu.cca <- readRDS("~/models/labelTransferCCA_union_hvg_F74_SCElist_20191113.RDS")
seu.liger <- readRDS("~/models/labelTransferLiger_union_hvg_F74_SCElist_20191113.RDS")
seu.conos <- readRDS("~/models/labelTransferConos_union_hvg_F74_SCElist_20191113.RDS")
integrate_features <- scan("~/intFeatures_union_hvg_2000_F74_SCElist_20191113.txt", what='')

Read 3029 items

int.list <- list(CCA=seu.cca, Liger=seu.liger, Conos=seu.conos)
# ## Make method color palette
# method.palette <- brewer_palette_4_values(names(int.list), "Set1")

Embeddings

Visualize label transfer on original ATAC data (embedded SnapATAC bins)

## Load original data
orig.ATAC <- readRDS("~/my_data/cellranger-atac110_count_30439_WSSS8038360_GRCh38-1_1_0.snapATAC.RDS")
sce.list <- readRDS("~/my_data/integrated_thymus/F74_SCElist_20191113.RDS")
orig.RNA <- sce.list$RNA
## Make SeuratObjects
atac.seu <- snapToSeurat(
    obj=orig.ATAC, 
    eigs.dims=1:20, 
    norm=TRUE,
    scale=TRUE
    )

Epoch: checking input parameters ... 
Non-unique features (rownames) present in the input matrix, making uniquePerforming log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Centering and scaling data matrix

  |                                                                              
  |                                                                        |   0%
  |                                                                              
  |==                                                                      |   3%
  |                                                                              
  |====                                                                    |   6%
  |                                                                              
  |======                                                                  |   9%
  |                                                                              
  |========                                                                |  12%
  |                                                                              
  |===========                                                             |  15%
  |                                                                              
  |=============                                                           |  18%
  |                                                                              
  |===============                                                         |  21%
  |                                                                              
  |=================                                                       |  24%
  |                                                                              
  |===================                                                     |  26%
  |                                                                              
  |=====================                                                   |  29%
  |                                                                              
  |=======================                                                 |  32%
  |                                                                              
  |=========================                                               |  35%
  |                                                                              
  |============================                                            |  38%
  |                                                                              
  |==============================                                          |  41%
  |                                                                              
  |================================                                        |  44%
  |                                                                              
  |==================================                                      |  47%
  |                                                                              
  |====================================                                    |  50%
  |                                                                              
  |======================================                                  |  53%
  |                                                                              
  |========================================                                |  56%
  |                                                                              
  |==========================================                              |  59%
  |                                                                              
  |============================================                            |  62%
  |                                                                              
  |===============================================                         |  65%
  |                                                                              
  |=================================================                       |  68%
  |                                                                              
  |===================================================                     |  71%
  |                                                                              
  |=====================================================                   |  74%
  |                                                                              
  |=======================================================                 |  76%
  |                                                                              
  |=========================================================               |  79%
  |                                                                              
  |===========================================================             |  82%
  |                                                                              
  |=============================================================           |  85%
  |                                                                              
  |================================================================        |  88%
  |                                                                              
  |==================================================================      |  91%
  |                                                                              
  |====================================================================    |  94%
  |                                                                              
  |======================================================================  |  97%
  |                                                                              
  |========================================================================| 100%

atac.seu <- RenameCells(atac.seu, new.names = orig.ATAC@metaData$barcode)
## Add cell type predictions
getPredictedLabels <- function(seu.int, int.name, id.col="predicted.id", score.col="score"){
  pred.df <- seu.int$ATAC@meta.data[,c(id.col, score.col), drop=F] 
  rownames(pred.df) <- str_remove(rownames(pred.df), "^ATAC_")
  colnames(pred.df) <- c(str_c("predicted.id", "_", int.name), str_c("score", "_", int.name))
  pred.df
  }
pred.cca <- getPredictedLabels(seu.cca, "CCA", score.col = "prediction.score.max")
pred.liger <- getPredictedLabels(seu.liger, "Liger")
pred.conos <- getPredictedLabels(seu.conos, "Conos")
if (all(rownames(pred.conos) == rownames(pred.cca)) & all(rownames(pred.conos) == rownames(pred.liger))) {
  atac.seu <- AddMetaData(atac.seu, metadata = cbind(pred.cca, pred.liger, pred.conos))
} else {
  stop("Non corresponding cell names")
}

plotly::ggplotly(pl)

geom_GeomTextRepel() has yet to be implemented in plotly.
  If you'd like to see this geom implemented,
  Please open an issue with your example code at
  https://github.com/ropensci/plotly/issues

Prediction score

Quantifies the uncertainty of the prediction. Calculated differently for every method, but used to define which cells are “unassigned”.

orig.composition <- orig.RNA$annotation
orig.frac <- table(orig.composition)/length(orig.composition)
orig.frac.df <- data.frame(orig.frac) %>%
  dplyr::rename(predicted.id=orig.composition, frac.label=Freq) %>%
  mutate(method="original.RNA")
score_cols <- str_subset(colnames(atac.seu@meta.data), 'score_')
label_cols <- str_subset(colnames(atac.seu@meta.data), 'predicted.id_')
pred.labels.df <- imap(list(CCA=pred.cca, Liger=pred.liger, Conos=pred.conos), ~ 
      rownames_to_column(.x, "cell") %>%
      rename_all(funs(str_remove(., str_c("_",.y)))) %>%
      mutate(method=.y)
    ) %>%
  purrr::reduce(bind_rows) %>%
  mutate(score=ifelse(is.na(score), 0, score))

funs() is soft deprecated as of dplyr 0.8.0
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once per session.binding character and factor vector, coercing into character vector

predict_score_hist <- 
  pred.labels.df %>%
  ggplot(aes(score, fill=method)) +
  geom_histogram(position="identity", alpha=0.8, bins=40) +
  facet_grid(method ~.) +
  scale_fill_brewer(palette="Set1") +
  xlab("Label prediction score") +
  theme_bw(base_size = 16) +
  theme(legend.position = "top")
cutoffs <- seq(0,1,0.05)
predict_score_cumedist <-
  pred.labels.df %>%
  group_by(method) %>%
  mutate(bins=cut(score, breaks = cutoffs)) %>%
  mutate(score=as.numeric(str_remove_all(as.character(bins), ".+,|]"))) %>%
  ggplot(aes(score, color=method)) +
  stat_ecdf(size=0.8, alpha=0.7) +
  scale_color_brewer(palette = "Set1") +
  ylab("Fraction of unassigned cells") +
  xlab("Prediction score cutoff") +
  theme_bw(base_size = 16) +
  xlim(0,1) +
  coord_fixed() +
  guides(color="none") 
ggpubr::ggarrange(predict_score_hist, predict_score_cumedist, common.legend = TRUE, widths = c(0.8, 1.2),
          labels=c("A", "B")) +
  ggsave(paste0(outdir, "prediction_score_distribution.png"), height = 6, width = 10)

Removed 63 rows containing non-finite values (stat_ecdf).

ggpubr::ggarrange(
  plotlist = list(
    FeaturePlot(atac.seu, reduction = "umap.snap", feature = "score_CCA"  , coord.fixed = TRUE) + ggtitle("CCA"),
    FeaturePlot(atac.seu, reduction = "umap.snap", feature = "score_Liger", coord.fixed = TRUE) + ggtitle("Liger"),
    FeaturePlot(atac.seu, reduction = "umap.snap", feature = "score_Conos", coord.fixed = TRUE) + ggtitle("Conos")
  ),
  common.legend = TRUE, ncol=3, nrow=1
) +
  ggsave(paste0(outdir, "prediction_score_umaps.png"), height = 7, width=14)

Cell type composition

Compare cell type fractions (w uncertainty)

pred.labels.df %>%
  group_by(method) %>%
  drop_na() %>%
  mutate(tot.cells=n()) %>%
  ungroup() %>%
  group_by(method, predicted.id) %>%
  summarise(tot.label = n(), tot.cells = max(tot.cells), mean.score=mean(score)) %>%
  mutate(frac.label=tot.label/tot.cells) %>%
  bind_rows(orig.frac.df) %>%
  mutate(orig.rank = orig.rank.df[predicted.id,]) %>%
  mutate(predicted.id=factor(predicted.id, levels=rownames(orig.rank.df)))%>%
  # select(method, predicted.id, frac.label) %>%
  # distinct() %>%
  ggplot(aes(predicted.id, frac.label, fill=mean.score, color=mean.score)) +
  geom_point(size=2) +
  geom_col(width=0.05) +
  coord_flip() +
  # geom_line(aes(group=method)) +
  facet_wrap(method~., nrow=1, ncol=4, scales="free_x") +
  scale_color_viridis_c() +
  scale_fill_viridis_c() +
  ylab("Fraction of cells") +
  theme_bw(base_size = 16) +
  ggsave(paste0(outdir, "cell_type_composition_bars.png"), width = 15, height = 7)

binding character and factor vector, coercing into character vector

Agreement with unsupervised clustering of ATAC data

Calculate which fractions of NNs in bin based graph of ATAC cells have the same annotation

k = 30
atac.seu <- FindNeighbors(atac.seu, assay = "ATAC", reduction = "SnapATAC", dims = 1:15, k.param = k)

Computing nearest neighbor graph
Computing SNN

atac.nn.list <- getNNlist(atac.seu)
score.CCA <- imap_dbl(atac.nn.list, ~ sum(pred.cca[.x,1] == pred.cca[.y,1])/k) %>% setNames(names(atac.nn.list))
score.Conos <- imap_dbl(atac.nn.list, ~ sum(pred.conos[.x,1] == pred.conos[.y,1])/k) %>% setNames(names(atac.nn.list))
score.Liger <- imap_dbl(atac.nn.list, ~ sum(pred.liger[.x,1] == pred.liger[.y,1])/k) %>% setNames(names(atac.nn.list))
knn_score_df <-
  as.data.frame(cbind(score.Conos, score.Liger, score.CCA)) %>%
  rownames_to_column("cell") %>%
  pivot_longer(cols=str_subset(colnames(.), "score"), names_to = "method", values_to = "KNN_score") %>%
  dplyr::mutate(KNN_score=ifelse(is.na(KNN_score), 0, KNN_score),
                method=str_remove(method, "score."))
quants = seq(0,1, by = 0.05)
AUECDF_knn_score <- knn_score_df %>%
  split(.$method) %>%
  map_dbl( ~ .x %>%
      arrange(KNN_score) %>% 
      {ecdf(.$KNN_score)(quants)} %>% AUC(quants,.)
    )
  
knn_score_df %>%
  mutate(AUC=AUECDF_knn_score[method]) %>%
  ggplot(aes(KNN_score, color=method, fill=method)) +
  stat_ecdf(size=1) +
  scale_color_brewer(palette = "Set1") +
  geom_text(data=. %>% group_by(method) %>% summarise(AUC=max(AUC)), 
            x=0.05, hjust=0,
            aes(label=glue("AUC = {round(AUC, 3)}"), y=c(0.90, 0.95, 1))) +
  theme_bw(base_size = 16) +
  ylab("ECDF") +
  ggsave(paste(outdir,"KNN_score_ecdf_unionHVG.png"), height = 4, width=6)

full_join(pred.labels.df, knn_score_df) %>%
  ggplot(aes(KNN_score, color=method)) +
  stat_ecdf() +
  facet_wrap("predicted.id") +
  scale_color_brewer(palette = "Set1") +
  coord_fixed()

Joining, by = c("cell", "method")

Accessibility of markers

Taking markers from Fig. S2 of JP’s manuscript

Reproducing Fig.2H on T-cell development

tcells.markers.df %>%
  full_join(t.cell.markers.df) %>%
  # filter(method=="CCA") %>%
  mutate(predicted.id=factor(predicted.id, levels=ordered.tcells)) %>%
  ggplot(aes( predicted.id, gene)) +
  facet_grid(cell.type.class~method, scales = "free_y", space="free") +
  geom_point(aes(size=frac.cells, color=mean.acc)) +
  scale_color_gradient(high="darkblue", low="white") +
  # scale_color_gradient2(midpoint = 0.5) +
  theme_bw(base_size = 16) +
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5),
        strip.text.y = element_text(angle=0)) 

Joining, by = "gene"
Column `gene` joining factor and character vector, coercing into character vector

Thoughts

• Conos scores a lot of cells with high confidence, but fails to assign cells to difficult clusters
• CCA resembles the composition of the RNA data better, but curious that the other methods identify way more